home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 November / macformat-018.iso / Utility Spectacular / Utilities / Calc / lib / poly.cal < prev    next >
Encoding:
Text File  |  1992-02-24  |  3.5 KB  |  211 lines  |  [TEXT/????]
  1. /*
  2.  * Copyright (c) 1992 David I. Bell
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * Calculate with polynomials of one variable.
  7.  */
  8.  
  9. obj poly {deg, coef};
  10.  
  11.  
  12. define pol()
  13. {
  14.     local x, d, i;
  15.  
  16.     d = param(0) - 1;
  17.     if (d < 0)
  18.         quit "No coefficients for pol";
  19.     if (d == 0)
  20.         return param(1);
  21.     obj poly x;
  22.     x.deg = d;
  23.     mat x.coef[d+1];
  24.     for (i = 0; i <= d; i++)
  25.         x.coef[d-i] = param(i+1);
  26.     return x;
  27. }
  28.  
  29.  
  30. define poly_print(a)
  31. {
  32.     local i, n;
  33.  
  34.     for (i = a.deg; i >= 0; i--) {
  35.         n = a.coef[i];
  36.         if (n == 0)
  37.             continue;
  38.         if (i == a.deg) {
  39.             if (isreal(n) && (n < 0)) {
  40.                 print "- " : ;
  41.                 n = abs(n);
  42.             }
  43.         } else {
  44.             if (!isreal(n) || (n > 0))
  45.                 print " + " : ;
  46.             else {
  47.                 print " - " : ;
  48.                 n = abs(n);
  49.             }
  50.         }
  51.         if ((n != 1) && (i > 0)) {
  52.             if (isreal(n))
  53.                 print n : "*" : ;
  54.             else
  55.                 print "(" : n : ")*" : ;
  56.         }
  57.         switch (i) {
  58.             case 0:
  59.                 if (isreal(n))
  60.                     print n : ;
  61.                 else
  62.                     print "(" : n : ")" : ;
  63.                 break;
  64.             case 1:
  65.                 print "X" : ;
  66.                 break;
  67.             default:
  68.                 print "X^" : i : ;
  69.         }
  70.     }
  71. }
  72.  
  73.  
  74. define poly_add(a, b)
  75. {
  76.     local x, d;
  77.  
  78.     if (isnum(b)) {
  79.         x = a;
  80.         x.coef[0] += b;
  81.         return x;
  82.     }
  83.     if (isnum(a)) {
  84.         x = b;
  85.         x.coef[0] += a;
  86.         return x;
  87.     }
  88.     if (a.deg == b.deg) {
  89.         d = a.deg;
  90.         while (a.coef[d] == -b.coef[d])
  91.             if (--d <= 0)
  92.                 return a.coef[0] + b.coef[0];
  93.     }
  94.     d = max(a.deg, b.deg);
  95.     obj poly x;
  96.     x.deg = d;
  97.     mat x.coef[d+1];
  98.     while (d >= 0) {
  99.         if (d > a.deg)
  100.             x.coef[d] = b.coef[d];
  101.         else if (d > b.deg)
  102.             x.coef[d] = a.coef[d];
  103.         else
  104.             x.coef[d] = a.coef[d] + b.coef[d];
  105.         d--;
  106.     }
  107.     return x;
  108. }
  109.  
  110.  
  111. define poly_neg(a)
  112. {
  113.     local x, i;
  114.  
  115.     x = a;
  116.     for (i = x.deg; i >= 0; i--)
  117.         x.coef[i] = -x.coef[i];
  118.     return x;
  119. }
  120.  
  121.  
  122. define poly_sub(a, b)
  123. {
  124.     return a + (-b);
  125. }
  126.  
  127.  
  128. define poly_mul(a, b)
  129. {
  130.     local x, i, j;
  131.  
  132.     if (isnum(b)) {
  133.         if (b == 0)
  134.             return 0;
  135.         if (b == 1)
  136.             return a;
  137.         if (b == -1)
  138.             return -a;
  139.         x = a;
  140.         for (i = x.deg; i >= 0; i--)
  141.             x.coef[i] *= b;
  142.         return x;
  143.     }
  144.     if (isnum(a)) {
  145.         if (a == 0)
  146.             return 0;
  147.         if (a == 1)
  148.             return a;
  149.         if (a == -1)
  150.             return -a;
  151.         x = b;
  152.         for (i = x.deg; i >= 0; i--)
  153.             x.coef[i] *= a;
  154.         return x;
  155.     }
  156.     obj poly x;
  157.     x.deg = a.deg + b.deg;
  158.     mat x.coef[x.deg+1];
  159.     for (i = a.deg; i >= 0; i--)
  160.         for (j = b.deg; j >= 0; j--)
  161.             x.coef[i+j] += a.coef[i] * b.coef[j];
  162.     return x;
  163. }
  164.  
  165.  
  166. define poly_div(a, b)
  167. {
  168.     local i, x;
  169.  
  170.     if (!isnum(b))
  171.         quit "Only division by numbers currently allowed";
  172.     if (b == 0)
  173.         quit "Division by zero";
  174.     if (b == 1)
  175.         return a;
  176.     if (b == -1)
  177.         return -a;
  178.     x = a;
  179.     for (i = x.deg; i >= 0; i--)
  180.         x.coef[i] /= b;
  181.     return x;
  182. }
  183.  
  184.  
  185. define ev(a, x)
  186. {
  187.     local i, r;
  188.  
  189.     obj poly r;
  190.     if (!istype(a, r))
  191.         quit "Evaluating non-polynomial";
  192.     i = a.deg;
  193.     r = a.coef[i];
  194.     while (--i >= 0)
  195.         r = r * x + a.coef[i];
  196.     return r;
  197. }
  198.  
  199. global lib_debug;
  200. if (!isnum(lib_debug) || lib_debug>0) print "obj poly {deg, coef} defined"
  201. if (!isnum(lib_debug) || lib_debug>0) print "pol() defined"
  202. if (!isnum(lib_debug) || lib_debug>0) print "poly_print(a) defined"
  203. if (!isnum(lib_debug) || lib_debug>0) print "poly_add(a, b) defined"
  204. if (!isnum(lib_debug) || lib_debug>0) print "poly_neg(a) defined"
  205. if (!isnum(lib_debug) || lib_debug>0) print "poly_sub(a, b) defined"
  206. if (!isnum(lib_debug) || lib_debug>0) print "poly_mul(a, b) defined"
  207. if (!isnum(lib_debug) || lib_debug>0) print "poly_div(a, b) defined"
  208. if (!isnum(lib_debug) || lib_debug>0) print "ev(a, x) defined"
  209. if (!isnum(lib_debug) || lib_debug>0) print "Use pol() to make polynomials (high coefficient first)"
  210. if (!isnum(lib_debug) || lib_debug>0) print "Use ev(a, x) to evaluate them"
  211.